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ABSTRACT 

Calculation of unsteady flow phenomena requires careful attention to the numer- 
ical treatment of the governing partial differential equations. The personal computer 
provides a convenient and useful tool for the development of meshes, algorithms, and 
boundary conditions needed to provide time accurate solution of these equations. The 
one-dimensional equation considered in this paper provides a suitable model for le 
study of wave propagation in the equations of transonic small disturbance potential 
flow. Numerical results for effects of mesh size, extent, and stretching, time step size, 
and choice of far-field boundary conditions are presented. Analysis of the discretized 
model problem supports these numerical results. Guidelines for suitable mesh and time 
step choices are given. 


1-INTRODUCTION 

At present, transonic small disturbance (TSD) potential theory provides the only 
practical method for the aeroelastic analysis of complete aircraft configurations at tran- 
sonic speeds (Batina et al. [1]). hi order to compute a time accurate flow field, careful 
attention must be given to the mesh design and boundary conditions Mesli fineness, 
extent, and stretching all affect the accuracy with which the wave-like disturbances 
inherent in the unsteady flow field can be described. 

In applications of an early TSD code (Whitlow (2]), problems associated with poor 
mesh design were encountered in calculating the unsteady lift aiid moment on an oscil- 
lating airfoil. The nature of the inaccuracies is illustrated in figure 1 (Seidel [3]). Shown 
are the real and imaginary parts of the lift curve slope for a flat-plate airfoil as func- 
tions of pitching frequency, compared with the exact linear theory results (Bland [4|) 
with which they should agree. Figure 1(a) shows results calculated using an expo- 
nentially stretched mesh which extends several hundred chord-lengths from the airfoil. 
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2-ANALYSIS AND RESULTS 

The transonic small disturbance equation in its simplest form may be written as 

(1 - M 2 - (7 + 1 )M 2 <j> £ )4> xx + 4 > zz = 2A i 2 <j> xt + M 2 tj) tt (!) 

floW - Elementar y solutions of this equation for Unear subsonic flow 
present circular waves radiating outward as their centers are converted downstream 
Waves propagating along a particular direction may be investigated by considering one' 
dimensional wave motion. In particular the difficulties encountered in the TSD results 
of figure 1 have been identified with wave reflections in the vertical direction 

2*1— One-dimensional wave equation 

The effects of mesh extent, spacing, and stretching on the solution of the TSD flow 
equation are clearly shown in the study of the wave equation 


<f>tt = c 2 <j> lz 


( 2 ) 


This equation is obtained horn equation ( 1 ) by omitting the ^-derivative terms and 

r? " Pl 8 y j. ‘ he r ° le °< of wave propagation speed c m 

the nondmiensional coordinates used. ” 

The usefulness of this approach in understanding problems encountered in TSD 
calculations is demonstrated in figure 2 . Numerical solutions are given for the response 
to a pulse m <j> (analogous to the pulse in angle of attack used to generate the TSD 
es s o gure ). each case the finite-difference potential y?( 0 ,<) at 2 = 0 is given 
as a function of time m the left part of the figure. The resulting frequency response 
unctions F(k) are plotted on the right. Figures la and 2 a use the same exponentially 
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stretched mesh in the z-direction. The quadratically stretched mesh of 25 chord-lengths 
extent proposed by Seidel [3] is used in figures lb and 2b-c. Figure 2c shows the effect of 
using an improper far-field boundary-condition. The results of figure 2 are remarkably 
similar to those shown in figures 1-2 of Edwards [6] for the TSD equation. In particular, 
the oscillations in the time history in figure 2a are known to originate from the 

large mesh stretching used and typically lead to oscillations in the frequency response 
function F{k) similar to those seen in this figure. Furthermore, the disturbance in the 
time history at t ~ 50 in figure 2c originates from a far-field boundary reflection and 
leads to the oscillations seen in F(k) at low values of k. 

The wave equation (2) has a general solution of the form 

<j>(z,t) = F(z - ct) + G(z + ct) (3) 

in which F and G represent, waves propagating in the positive (z — ct constant) and 
negative (z + ct constant) z-directions with speed c. 

To simulate the vertical component of the fluid flow problem, we may solve this 
problem in the quarter plane z,t > 0 subject to two initial conditions at t = 0 and 
boundary conditions at z = 0 and oo. To solve our model problem numerically, we must 
discretize the continuous variables z, t and replace the partial derivatives in equation (2) 
by finite-difference analogues. Finally, the boundary condition at 2 = oo must be 
replaced by a condition at some finite value z = L. 

Before proceeding with the details of the numerical approach, several illustrations 
of the nature of the difficulties with numerical solutions of an initial value problem 
will be given. Calculations for an initial wave packet at t = 0 are shown in figures 3. 
Exponentially stretched meshes of varying fineness (illustrated beneath each figure) are 
used and the solutions are given at every other time step using a small step size, hi 
figure 3(a) the wave propagating outwards encounters a reflecting boundary at the right 
and a non-reflecting boundary at the left. The initial pulse divides into two pulses 
which propagate outward with little distortion until the boundary is encountered, hi 
figure 3(b) the same 97-point mesh is used with non-reflecting boundary conditions at 
both sides. Figure 3(c) is for a coarser mesh which allows significant reflection from the 
boundaries to occur, hi figure 3(d) the mesh has become so coarse that essentially all 
of the wave is reflected from within the mesh, even before reacliing the boundaries. For 
this case, at the boundaries there are only about two mesh points per wave length. Most 
of the wave energy is trapped within the mesh. In an aerodynamic problem, the wing 
would remain near z = 0 and would encounter these non-physical mesh reflections. 
The initial pulse contains components of all frequencies and, as we shall see, the liiglier 
frequencies propagate more slowly and therefore lag behind the main (low frequency) 
wave. 

The initial conditions for the remainder of the paper will be 

<£(z, 0) = 0 and ^ f (z,0) = 0 (4) 

As the boundary condition at z = 0 we use 


«M<M) = /(<) 


( 5 ) 


This condition is analogous to imposing a condition on the vertical velocity (downwash) 
at the oscillating airfoil surface in the aerodynamics problem. Finally, we must impose 
a boundary condition at x = L which simulates the outgoing wave condition at infinity 
m the original problem. The proper choice for this wave propagating upward is 

t) + c<j> z (L, t) = 0 (g) 


2.2-Finite-difFerence formulation 

The one- dimensional model equation (2) is discretized in a manner identical to that 
used for the vertical * and time t coordinates of the TSD equation in the program 
described by Batina [7]. The second-order accurate backward difference approximation 
for the t-derivative at the new time level n -f 1 is 


d n <f = (2ip n+1 - 5 <p n + 4 v? n - 1 - <p n ~ 2 )/At 2 (7) 

winch involves four time levels. For reference we also state the first-order accurate form 

dtt? = (v> n+1 - 2<p n + <f n ~ l )/At 2 (8) 

The second-order accurate central difference approximation for the 2 -derivative is 


d zzV > = - ( , ii±LZh 

1 ~ z j~ l ' ~ z j z j ~~ z j-i ) 

in which this equation applies at the internal mesh points 0 - j < J 

The boundary condition in equation (5) at 2 = 0 may be incorporated 
the differential equation (2) at 2 0 using 


(9) 


by solving 



in which the zeroth mesh point need be only near 2=0. 

The far-field condition in equation ( 6 ) is derived by applying the boundary condition 
a * j = J - 1/2. The second-order accurate formulas needed are 


~J - ZJ - 1 


and 

dt<f = (3^ n+1 - 4y>" + <p n ~ 1 )/2At 
and the first-order accurate formula is 


( 11 ) 

( 12 ) 


dt<P = (v ?nfl - (13) 

The ^-derivatives in equations (12) and (13) are evaluated at j = J - 1/2 by applying 
the operations on the right hand sides to the average (yj-i + <pj)/2. 
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2.3-Model equation 

There is no loss of generality in normalizing the time t such that the wave speed 
becomes c = 1 . This choice slightly simplifies the analysis. 


Continuous Problem To be specific, in the quarter-plane t,z > 0, solve 


<}>tt — fizz 

for t, z > 0 

(14a) 

<f>{z, 0) = o) = 0 

for initial conditions 

(14b) 

<M<M) = f(t) = -u>sinu )t 

for downwasli b.c. 

(14c) 

<f> t + <f> z = 0 as z — * oo 

for outgoing wave 

(14d) 

The exact solution is 



<f>{z,t) = 1 — cosa>(< — z) 

for z < t 

(15) 

and 0 otherwise. This solution represents an undistorted wave propagating upward with 
unit speed c = 1 . We remark that this solution satisfies the outgoing wave condition of 


equation (14d) everywhere. 


Semi-Discrete Problem The semi-discrete analysis (‘method of fines’) considers the 
effect of the spatial discretization only, while treating the time as continuous. The mesh 
characteristics are thereby isolated from the questions associated with time step size 
(i.e., Courant number v) such as stability. 

For a uniform mesh 


Zj = hj for j = 0, 1, . . . , J (16) 

with constant spacing h = A z, central differencing provides 0(h 2 ) accuracy. Letting 
ip = dip/dt, the semi-discrete version of equation (14) is 

if 0 = (2/h 2 ){(pi - y?o) + (2u>/h) smart (17a) 

= (l/h 2 )(<fj-i - 2ipj + (17i>) 

<fj-i + = (2/h)(<pj-.i - <pj) (17c) 

Vichnevetsky and Bowles [8] and Vichnevetsky [5, 9], provide an extensive treat- 

ment of the problems encountered in using finite-difference methods to solve equations 
governing wave-like phenomena. They concentrate on the first order advection equa- 
tion (a model for the Euler equations) rather than the wave equation treated herein. 
Trefetlien [10] gives the analysis and some numerical results for both the one- and two- 
dimensional equations. For the analysis, we assume that the solution is harmonic in 
time and space and determine the conditions implied by this assumption. For more com- 
plicated wave-forms the ideas of Fourier analysis apply. This approach does not capture 
the transient effects present in the initial value problem, but looks at the steady-state 
behavior. 

Substitution of the harmonic trial solution 

^j(t) = e i(w ‘" < ** ) 


(18) 



into e<lualiou (17b) yields the dispersion relation. Dispersion refer, to the fact that 
wares with different frequencies will propagate at different speeds. We obtain 

. £h wh 

sm T = ± T (19) 

in which we consider the temporal frequency u as given and the spatial ware number f 
*s unl[now “- Choosing the + sign for the upward propagating ware yields ' 


> 2 , ivh 

4=7 arc sin — 
h 2 


w ( 1 + 


1 / uh \ 2 \ 
0 \~2 ) ) 


> ui 


( 20 ) 


Note that as u>h 0, we recover, as expected, the exact solution for which ( = u>. 

If we restrict attention to the outgoing wave solution (semi-infinite mesh extent ) the 
+ sign m equation (19) is selected and imposing the boundary condition equation (17a) 
hen yields the exact solution for the semi-discrete problem 


<M<) = i - 


cos(u>< - £zj) 


y/l^ {u>h/ 2) 2 ‘ 

wliich is also an exact solution of the ‘continuous’ problem 


( 21 ) 


Vtt = c*V« 


in wliich the phase velocity 


( 22 ) 


r = 




£ arcsinoi/i/2 


6 


( 


2 ) 


(23) 


replaces the original wave speed c = 1, provided that we identify the discrete variable z, 
with the continuous variable z. 

1 J he V: SUl l S proV j de a great deal of insight, into the behavior of the numerical 
solution. Note that it depends only on c*, a function of so that increasinTehW 

hi 1 f V Same C , EqUati ° n ( 19 ) M/2| < 1 (cut-off frequency) 

m order for harinmuc solutions to exist, although other types of behavior may appear 

( Viclmevetsky [5, 9)). The wave propagates with the phase velocity c* and wave length 


A* = 2n/( = C *X < X 


(24) 


where w - < £ < n/h so that the numerical wave propagates more slowly than the exact 
wave and has a shorter length. The amplitude is too large by the factor 


A = 


1 

y/1 - M/2) 2 




(25) 


The decreasing propagation speed c* 
of the strange behavior of the solution 


with increasing mesh spacing h lies at the root 
on the stretched meshes illustrated in figure 3. 
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The analysis shows that these effects are present in the semi-discrete equation and do 
not depend on the time step used. 

The decrease in propagation speed (equation (23)) and increase in amplitude (equa- 
tion (25)) for a harmonic wave are clearly seen in figure 4. In figure 4(a) the solution 
<p(z,t) on a uniform 49-point mesh is shown at t = 16. The numerical wave is prop- 
agating more slowly across the mesh than the exact wave. In figure 4(b) the solution 
time history y>(0,<) is shown. For the parameters used here, u> = tt, h = 1/3, we 
expect c* = 0.95 and A = 1.17. The numerical results agree with these values. For 
this example there are six mesh points per spatial wave length, which is a reasonable 
minimum number to capture the propagating signal. 

The analysis considered to this point correctly describes the behavior of the numeri- 
cal solution of the semi-discrete equation on a semi-infinite mesh as t -> oo. Considering 
now the effect of limited mesh extent, we will show that the second-order accurate non- 
reflecting far-field boundary condition equation (17c) produces a partial reflection of 
the outward propagating wave. For an incident wave, equation (21) with amplitude A 

<p l = 1 — Acos(u;< — £z) (26) 


the reflected wave is 

<p r = —A cos(u>f + £(z — 2 L)) (27) 

A ~f" 1 

where L = (zj_i + zj)/2 is the numerical boundary location. This reflected wave has 
an amplitude reduced from that of the incident wave by the reflection ratio 


A - 1 1 / uh \ 2 

P= A + 1 ~~ 4 V 2 ) 


(28) 


and is propagating in the negative 2 -direction from an image point located two mesh 
extents from the origin. It is difficult to observe such reflections for the harmonic 
boundary condition under discussion here. However, they will be seen in the results for 
the pulse boundary condition given below. 

The semi-discrete harmonic analysis may also be used to show that changes in 
mesh spacing produce internal (as opposed to boundary) reflections. This analysis lias 
application to the results observed on stretched meshes. Consider a mesli with a single 
change in spacing at : = 0. Let the spacing be h for z < 0 and H for z > 0 
with corresponding wave numbers £/, and £h- We wish to determine the reflection 
that occurs when the upward propagating wave of unit- amplitude passes through the 
discontinuity. We remark that the change in i at the interface leads to a change in the 
phase velocity c* there. Away from the interface the solutions will be 

Vj (t) = + pe' (u ' ,+0, " j) (29) 




for < 0 and 


( 30 ) 
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for zj > 0. Here p is the reflection ratio and r is the transmission ratio. These 
parameters are determined by applying equation (9) at the interface point Zj yielding 


Vi = 


H + h 


( Vj+i ~ 

V H 


Vj _ Vi zJLLzl 

h 


and by enforcing continuity to give a second equation 

T = l + P (32) 

The solution is a little complicated and it suffices to show the asymptotic results for 
small mesh spacmgs and/or frequency 

?-(j)V-fc 2 ) (32) 


r 1 + 


(i)V 


Note that as the wave passes through an interface with increased mesh spacing (H > 

k >' of the transmitted wave is increased. Tliis result suggests that the 

amplitude would continue to grow as a wave propagates through a stretched mesh like 
those normally used in solving exterior flow problems. 

These reflections are illustrated in figure 5 for a pulse in ^(0,/). The pulse is 
centered at t c % 2 and the mesh has spacings h = 3/8 for 0 < z < 6 and H = 3/4 
or 6 * < 12 ( bi-uniform mesh’). In figure 5(a) the exact and numerical waves are 

shown at / 6. The solutions at the boundary ^(0,<) are shown in figure 5(b). The 

bold tic-marks indicate the times for the center of the pulse and its reflection from the 
f.ir-held to occur. Three different effects may be noted: (1) the return of the boundary 
reflection at t 2r 26; (2) the reflection generated by the mesh discontinuity at t % 14- 
(3) the oscillations near t * 3 wliicli are caused by too large a spacing to accurately 
define the pulse boundary condition imposed at r = 0. hi figure 5(c) the effects of these 
inaccuracies on the frequency response function are shown. Prom the view point of the 
aeroelasticiaii, the results might be satisfactory since little error occurs for k < 0.3. 

Additional insight into the behavior of the solution of the semi-discrete equa- 
lon ( 31) on stretched meshes may be obtained by examining the modified equation 
( n erson et al. [11]). Tins equation is obtained by Taylor series expansion of the 
potentials. For the uniform mesh the modified equation contains only even derivative 
erms and has constant coefficients. For the nonuniform mesh, derivatives of all orders 
are present and, in general, the coefficients are functions of r. A unique feature of the 
quadrat, cally stretched mesh proposed by Seidel [3] is that the coefficient of the third 
derivative term in the modified equation is a small constant. 

Fully Discrete Problem For simplicity, the analysis of the fully discrete equations 
will be given only for the case of the uniform mesh. As in Batina [7], the time differencing 
is written in implicit form. Both the first- and second-order time accurate equations 
w,U be displayed, then the dispersion relations will be given and the implications for 
the numerical solutions will be addressed. 


(34a) 

(346) 

(34c) 


The equations of 0(A<), obtained from equations (8), (13), and (17) are 

(i + 2 i/ 2 k +i - 2 I /vr M = 2 - rr 1 - 2 hv 2 m 
f + <1 + 2^ 2 = 2 ^ - v-r 1 

(1 - 2v)<p n j±\ + (1 + 2 v)<p n j +1 = + ¥>" 


in which equation (34b) applies for 0 < j < J. 

The equations of 0(Af 2 ), obtained from equations (7), (12), and (17) are 

2(1 + t/ 2 Vo + 1 - 2i/V? +1 = 5v>0 - vr 1 + ¥>r 2 - *hv 2 f(t) (35 a) 
+ 2(1 + ^ 2 )^" +1 “ ^ 2 Vj+i = &Pj “ **7 + VT 2 ( 356 ) 

(3 - 4 vW£\ + (3 + 4i /)<^ +1 = 4(<^5_ a + ¥>j) - (V’jll + ^r 1 ) ( 35c ) 


in wliicli equation (35b) applies for 0 < j < J. 

The tridiagonal systems in equations (34) and (35), which include the boundary 
conditions, may be solved using the Thomas algorithm (Anderson et al. [11], p.128). 
These equations have been used for all of the numerical results presented herein. 

The dispersion relations are obtained by substituting the harmonic trial solution 

= e i(u ,t n -(zj) (36) 

with Zj = jAz and t n = nAt into the field equations and dividing by <?” • The result 
for the O(At) equation (34b) is 


1 . £Az 

Az Sm 2 “ A# 6 


1 A 4 /O A/ 

— tu>At/2 sin 


It suffices to give the asymptotic solution of this equation which is 
£ ~ u> (l + ^(Az 2 - 4Af2 ) ~ l '^ A# + °(^ 3A3 )) 


(37) 


(38) 


The result for the 0(Af 2 ) equation (35b) is 


1 . £A z 

—— sin — - — 

A z 2 


J_ -io,At/2 sin ^ y 2 _ e — iu/At 

At 2 


(39) 


for which 


£ - u. (l + A .- 3 + llAf 2 ) - A( J + 0 (id 4 A 4 )) 


(40) 


In both equations (38) and (40) the negative imaginary parts lead to an expo- 
nential decay of the solution with increasing z. hi particular the amplitudes decay 
as e VAc/2 for the 0 (At) equation and as e -^-At 3 z/i6 for the 0(A# 2 ) equation. 

Thus, in contrast with the result for the semi-discrete analysis, which predicts that the 
propagating harmonic wave will have its amplitude increased (see equation (25)), the 
effect of the time discretization is to produce both a spreading out of the frequencies and 
a decrease in amplitude with increasing distance z from the source of the disturbance 
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at the origin. Also, we see that the first-order method is much more dissipative than 
the second-order one. This greater dissipation makes the first-order method much more 
forgiving of poor mesh design and choice of far-field boundary conditions, at least in 
applications to problems for which only the near-field solution is of interest. 

The wave propagation speeds are computed from the real parts of the dispersion 
relations (equations (38) and (40)). For the two cases under consideration here, we see 
that the phase velocities will be approximately 

c* = Re{a>/£} - 1 - ^-(Az 2 - 4At 2 ) (41) 

for the O(Af) equation (38) and 

c* = Re{a>/£} - 1 - ^-(Az 2 + 11A< 2 ) (42) 

for the 0(A* 2 ) equation (40). Note that in equation (41) the error will vanish for 
v = Af/Az = 1/2, while no such possibility holds for the 0(At 2 ) equation. 

hi order to treat, a nonuniform mesh using the uniform mesh results we may apply 
the conclusions just reached in a local sense ( Viclmevetsky [9]). Then, for example, as a 
wave propagates into a stretching mesh we would expect both the dispersion (separation 
of frequencies) and dissipation (amplitude decay) to increase as the local Az increases. 
One can think of the harmonic solution given above as remaining approximately valid 
by simply treating Az as a function of z. An alternative approach is to consider 
perturbation solutions of the nonconstant coefficient modified equation. 

The effects of time step size and order of accuracy are illustrated in the final four 
figures. In order to reduce the appearance of far-field boundary reflections, a mesh is 
used which stretches quadraticafly from both boundaries. The mesh contains 25 points 
and extends to z = 8. The mesli spacing varies from Az = 0.03 at the ends to 0.64 in 
the center. The time step sizes are A t = 1/8 and 1/2. The calculations are carried to 
a total time of about 64 to allow all transients and reflections to die out. 

Figures 6-7 give results using the second-order accurate equation (35). For the 
larger time step in figure 6, the solution spreads out as it propagates across the mesh 
(figure 6(a)). No reflections are apparent in the time history shown in figure 6(b). 
Although the frequency response results (figure 6(c)) are smooth, large errors are present 
for k 0.5 which would go un-noted were the exact results not available for comparison. 
For the smaller time step results of figure 7, an unanticipated difficulty appears in the 
many high frequency oscillations in the time history of figure 7(b). These oscillations 
are also seen in the trailing part of the spatial wave in figure 7(a). hi a TSD calculation, 
the appearance of oscillations like these in a solution might lead to nonlinear instability. 
However, the frequency response results of figure 7(c) are quite good for k < 1. 

Figures 8-9 give results using the first-order accurate equation (34). Here the 
results are quite smooth for both time step sizes and no inaccuracies are apparent in 
the frequency response results of either figure 8(c) or 9(c), although large errors occur 
for the larger time step (figure 8). Examination of figure 8(b) shows that the numerical 
result predicts that the pulse center occurs too early in time. An arbitrary shift of 
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A</2 in the application of the forcing function f(t) in equation (34a) corrects this 
discrepancy very nicely. 

In choosing a vertical mesh for use in solving the TSD equation, one would not 
in general be able to ‘waste’ mesh points at the far-field boundaries by employing the 
very small spacing used in these examples, especially in three- dimensions. Nevertheless, 
the spacing must not be allowed to grow too large in the mid-field and should reduce 
somewhat at the far-field boundaries in order to avoid the appearance of non-physical 
oscillations in the calculation. 

3- CONCLUSIONS 

The personal computer provides a convenient and useful tool for the study of the 
techniques of mesh and algorithm design needed to provide time accurate solutions of 
the unsteady flow equations. The one- dimensioned wave equation considered in this 
paper provides a suitable model equation for studying wave propagation in the vertical 
direction in transonic small disturbance potential flow. The following conclusions may 
be drawn from this study: 

(1) The appearance of non-pliysical oscillations in the solution may result from either 
(a) inappropriate use of reflecting boundary conditions, (b) too large a mesh spac- 
ing at the boundary to correctly model the boundary condition, or (c) too much 
stretching of the mesh, leading to internal reflections. 

(2) In order to accurately convect wave like disturbances through the finite-difference 
mesh, there must be at least six mesh points per spatial wave length. 

(3) First-order time accurate methods are much more forgiving of poor mesh design 
than are second-order methods. The higher dissipation present with the first-order 
methods reduces the appearance of unwanted reflections, but at the price of a much 
poorer simulation of the overall flow field. 

(4) Improved finite- difference meshes may be designed using the one- dimensional model 
problem as a guide. 

4- NOMENCLATURE 

A amplitude of wave 

c wave speed in exact solution 

c* phase velocity in discrete solution 

F(k) frequency response function 

f(t) boundary value of <f> z 

J index of last mesh point 

j index for discrete z 

k reduced frequency 

L far-field boundary for z 

M Mach number 

n index for discrete t 

t time 

x streamwise coordinate 

z vertical coordinate 
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A 
A* 

v 

i 

p 

T 

<t> 

<P 
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Figtire 1.- Unsteady forces calculated using TSD theory at M = 0.85. 



Figure 2.- Time history and frequency response solutions of the one- dimensional model 
problem using the meshes of the TSD solution in figure 1. 
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Figure 7.- Comparison of exact and second-order accurate numerical solutions with 
At 1/8. 
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